Stata：特征值、特征向量与主成分分析-pca
👇 连享会 · 推文导航 | www.lianxh.cn
🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享 💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归 ⛳ 专题：SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching 🔨 工具：工具软件 | Markdown | Python-R-Stata 🎧 课程：公开课-直播 | 计量专题 | 关于连享会
连享会 · 文本分析 | 爬虫 | 机器学习
作者：姚和平 (中山大学)
邮箱：yaohp@mail2.sysu.edu.cn
目录
1. 猫头鹰群落演化问题
2. 特征值与特征向量
2.1 相关定义
2.2 应用到动力系统
2.3 QR 算法
3. 主成分分析 (PCA)
3.1 基本原理
3.2 Stata 实操
4. 参考文献
5. 相关推文
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
特征值和特征向量是线性代数里十分重要的概念，在经济学、统计学、物理、化学等学科中有着广泛的应用。这篇文章将为读者细致地讲解这两个概念，并介绍其重要应用—PCA (Principal Component Analysis)。作为引入，我们先来看一个现实案例。
1. 猫头鹰群落演化问题
1990 年，环境保护学家与木材行业在砍伐太平洋西北部原始森林，是否会造成斑点猫头鹰灭绝的问题上争论不休。数学家们也因此开始了对斑点猫头鹰种群的动力学研究，R.Lamberson (2005) 及其同事建立了以每年种群数量为区间的演化模型，时间为 。猫头鹰的生命周期自然分为三个阶段：幼年期 (1 岁以前)、半成年期 (1-2 岁)、成年期 (2 岁以后)，第 年的种群量可以用向量 表示。
根据统计数据，新的幼年猫头鹰在 年的数量是成年猫头鹰在 年数量的 0.33 倍，每年 18% 的幼年猫头鹰得以生存进入半成年期，71% 的半成年猫头鹰和 94% 的成年猫头鹰生存下来被计为成年猫头鹰。
我们用数学公式来表达这一模型：
将上面式子改写成矩阵形式，就将该问题转换成一个形式为 的方程。要看猫头鹰会不会灭绝，只需要计算出 即可。
上式表达除了可以用来表示猫头鹰群落的演化过程，还可以表示类似城市人口流动、传染病 SEIR 模型等，在日常生活中应用十分广泛。但是 毕竟是一个矩阵的乘方运算，还要求极限，看上去非常的复杂和不直观，因此我们引入了特征值和特征向量的概念，它们能使上述计算变得十分简单。
2. 特征值与特征向量
2.1 相关定义
学过线性代数的读者都知道，矩阵乘法 的本质就是对向量 进行线性变换，但我们却很难想象到结果向量 与起始向量 有什么联系。它们并不像我们以前学的代数乘法那么形象，就是一个简单的倍数关系，原始量经过放大或缩小就可以得到结果。但有没有可能把复杂抽象向量乘法转化为普通的代数乘法呢？请看下面的例子 (Lay，2016)。
设 ，，。下图是 和 乘以 后的图像。可以看出 与 不同，正好就是 ， 就像普通乘法一样只是把 拉长了。
因此我们可以看到，对于一个矩阵 ，确实存在向量 使得向量乘法的效果与标量乘法等价，即 ，这样就得到了我们开始想要的结果。为了方便起见，人们把能使 的这个 称作矩阵 的特征值，向量 称为特征向量。
聪明的读者可能会想到，通过定义条件 并不太容易找到特征值和相应的特征向量，而且特征值和特征向量常常并不止一对。那我们是否有一个简单的方法来找到一个矩阵的全部特征值和它所对应的特征向量呢？
把 做一个小变形，变为 ，那么该式子的非零解 (也称为非平凡解) 就是我们想求得的特征向量。根据矩阵可逆定理，当且仅当 ， 存在非零解，于是我们可以得到推论：
数 是 矩阵 的特征值的充要条件是， 是特征方程 的根。
因此我们可以通过解特征方程 来求出一个矩阵的特征值，进而求出每个特征值对应的特征向量。以下将通过一个小例子来演示特征值求解过程 (Lay，2016)。
求 的特征值：
因式分解后得到的上式也被称为特征方程，因此 的特征值为 3 或者 -7，分别代回 就可以得出分别对应的特征向量是 (3, 1) 和 (1, -3)。
特征值和特征向量还有一个非常有趣的性质：
如果我们能够将矩阵 拆解成 ，而且矩阵 是一个对角矩阵 (除对角线元素外的元素均为零)，那么矩阵 中的列向量就是 的特征向量， 对角线元素就是 的特征值，而且各特征向量与特征值在列位置上相互对应。
这个过程也叫矩阵的对角化，此时矩阵 与 具有相同的特征值，我们称他们为相似矩阵。这个性质在之后讲解主成分分析 (PCA)中有重要作用。
2.2 应用到动力系统
我们回到斑点猫头鹰的演化问题，特征值与特征向量是解决动力系统的关键。你可能会问，这种将矩阵乘法变为标量乘法的方法对于任意一个特定方阵来说，只有少数几个特征向量有效，比如在斑点猫头鹰的问题里，种群量 并不一定是转移矩阵 的特征向量。答案是，我们可以将任意一个非特征向量 用矩阵 的特征向量的线性组合来表示 (只要那些特征向量线性无关)。
以斑点猫头鹰为例，这里 都是转移矩阵
的特征向量， 是猫头鹰种群 时的种群量。那么就有：其中 和 都是标量，依次迭代后一般有：
这样我们就把一个复杂的矩阵乘法问题转换为一个简单的标量乘法问题，但这样会对运算产生什么样的影响呢？具体到斑点猫头鹰问题，我们可以计算出转移矩阵 的特征值分别为 ，，。因为 ，所以三个特征值的绝对值都小于 1。因此我们可以将 时的猫头鹰种群量用下式表示出来：
由于所有的特征值都小于 ，所以当 趋向正无穷时，右边的每一项都趋于零向量。因此，实序列 也趋于零向量。很不幸，该模型最终预测斑点猫头鹰最终将全部灭亡 (Lamberson，1992；Lay，2016)。
2.3 QR 算法
对于求一个 矩阵特征值的问题，虽然上面我们介绍因式分解特征方程的算法。但在实际中，当 时，除非是精心选择的矩阵，否则对其进行因式分解并不容易，而且对于 的一般 矩阵，没有公式或有限算法来求解其特征方程。因此，最佳的求特征值的算法会完全避开特征多项式。
事实上，Matlab 的算法就是先计算出特征值 ，然后通过展开乘积 来得到 的特征多项式。那这种不需要解特征方程的方法是什么呢？
QR 算法的核心思路是由矩阵 出发产生一个矩阵序列，这些矩阵都相似于矩阵 ，而且一部分矩阵接近于上三角矩阵。线性代数中两条关于特征值的简单定理：
假如 和 是 矩阵, 如果存在可逆矩阵 ，使得 ，则称 相似于 ，相似的矩阵有相同的特征值； 上三角矩阵的对角线上的元素就是该矩阵的特征值。
根据上面两条定理，我们可以很容易看出， 产生的矩阵序列中的近上三角矩阵的对角线元素就是 的特征值 (当然这只是一个估算的结果，因为有时无法产生真正的上三角矩阵)。那么我们应该怎么产生这些相似的上三角矩阵呢？
QR 分解是指将矩阵 分解成两个矩阵 和 的乘积，其中 是正交矩阵，满足， 是一个上三角矩阵。Stata 中的 qrd() 函数可以完成这一过程。
// --QR分解--
mata //进入 Mata 环境
A = ( 1, 2, 3 \ 4, 5, 6 \ 7, 8, 9 )
qrd(A, Q, R)
st_matrix("Q",Q) //将 Meta 里的矩阵转换到 Stata 环境中
st_matrix("R",R) //同上
end //退出 Mata 环境
mat list Q
mat list R
当矩阵 拆解成两个矩阵 后，产生的新矩阵就是 ，紧接着我们又对 进行 分解，重复之前的步骤，这样每次产生的新矩阵 就是产生的矩阵序列。证明 与 相似是很简单的：
随着新矩阵 的不断产生，新矩阵会不断收敛成为一个上三角矩阵，这就是我们最终想要的结果 (因为上三角矩阵的对角线元素就是矩阵 的特征值)。将这个过程总结成算法是非常简洁的。
到目前为止，我们已经讲解清楚了特征值和特征向量，那么接下来让我们来看一下主成分分析 (Principal Component Analysis) ，它是特征值和特征向量的一个主要应用方向。
3. 主成分分析 (PCA)
在大数据时代，我们面对的数据往往更多的特征。过多的变量会使程序变得更慢，也会让我们更难以发现数据之间的关系。幸运的是，在真实的世界里，我们通常可以在不损失太多信息的前提下，把一些彼此关联十分紧密的变量进行合并和去掉一些无关变量。
例如，下面左图中的数据点看似在三维空间中但实际上都落在一个二维平面上，这是因为描述数据点的某一个变量与其他两个变量完全线性相关，因此只需要调整一下坐标系就可以去掉一个变量。右图中数字周围的空白部分不提供信息，同样可以去掉。通过对原数据进行一些调整来减少变量个数的方法就叫做降维，而 PCA 是迄今为止使用最广泛降维算法之一。
接下来，我们将介绍 PCA 的基本原理和 Stata 实操的内容。
3.1 基本原理
从上面的介绍中，我们可以看出降维算法主要有两个功能：去相关和去冗余。通俗来说，去相关是指尽可能减少变量间的相关性，去冗余是指去除那些方差很小的无关变量。PCA 算法可以完美地做到上述的功能，让我们来看一下它是怎么实现的。
假设现有数据集 ， 表示变量个数， 表示样本数量。那么有什么数学工具可以同时表示数据集中变量的相关性和方差呢？协方差矩阵！其定义如下，因此协方差矩阵对角线衡量变量间方差，其余元素衡量变量间的协方差。
协方差矩阵的求解公式如下：
为达到上述两个目标，我们希望对原数据集进行一个简单变换，使变换后矩阵的协方差矩阵除对角线外元素都等于 0，对角线上元素表示原来变量的方差。很自然我们会想到上文提到的对角化，因为对角矩阵恰好满足上述性质，因此将对原数据集的协方差矩阵进行对角化，得到其特征值和特征向量。
这时我们会发现矩阵 的协方差矩阵恰好为对角矩阵 ：
因此，用原数据集 的协方差矩阵的特征向量矩阵 对原矩阵进行变换，可以使变换后的数据集 的协方差矩阵等于原来的对角矩阵，即变换后的数据集 有如下性质：
变量间彼此不相关； 原数据集 协方差矩阵的特征值衡量了 各变量的方差。
我们将 中的特征向量称为主成分，对原数据集 右乘 相当于对原数据集变换了坐标轴，原来的坐标轴是各原始变量，新的坐标轴是各主成分。这种变换实现了这样的目的：在保证不损失信息的同时，使变换后各新变量不相关，原数据集 协方差矩阵的特征值大小表示各主成分解释新数据的方差大小。因此新数据集各主成分对应特征值的大小表示了该主成分的重要性，之后可以借此来对主成分进行筛选，达到去冗余的目的。
接下来 我们将进入降维的步骤，去除那些相对不重要的主成分。当去除一些主成分 (坐标轴) 时，我们相当于将数据从高维空间投影到一个低维空间。下图是一个简单的例子，现在将数据从二维平面投影到一条直线上，我们显然会选择投影到 上。因为这样保留了最大的方差 (信息量)，在主成分分析中我们也会做同样的选择。
因此在 PCA 算法中，我们只会保留那些投影后保留了原数据最大方差的主成分。而上式告诉我们，特征值越大，主成分方差越大。于是我们按照特征值从大到小的顺序将对应的特征向量从左到右排列，只保留变换矩阵 前 列，记为 ( 的大小可以自由选择，下面的实例讲解中会提供一些选择参考)。最终我们对上面的变换 修正为：
就是我们最终希望得到的降维后的数据矩阵。另一方面，上图也揭示了特征值、特征向量与主成分分析的核心联系。我们回到特征值的定义：
换一种理解方式，现在 是原数据集， 看成是将数据集 投影到向量上， 就衡量了数据集投影到某个特征向量 (主成分) 上的方差大小 (在上图里， 就是数据集投影到主成分后的长度)。我们已经知道 PCA 算法的本质就是寻找主成分，然后将数据集投影到主成分上的过程，而在 PCA 中主成分就是特征向量。
因此 PCA 算法的核心逻辑和特征值的定义式 有着惊人的吻合，特征值 的大小刚好就衡量了对应主成分解释原数据的方差，而这也正是理解特征值、特征向量与主成分分析之间关系的核心，这也是特征值的重要意义之一——它衡量了数据在特征向量方向上所包含的信息量！
值得注意的是，虽然我们上面用协方差矩阵进行推导，但在实际中我们常常用相关系数矩阵进行代替，两者作用相似，但前者比起后者更容易受到量纲的干扰。另一方面，并不是所有的数据都适合使用 PCA 算法。PCA 算法对数据变量之间的相关关系有限制，在使用前应先对数据进行检验。
3.2 Stata 实操
接下来，我们使用书中的一个例子来进行演示 (Jackson，2005)。该数据集记录了一百名成年男性左右耳的听力情况，测量结果是左耳和右耳在四个不同频率上的最小可辨识强度，比如变量 lft1000 指的是左耳在 1000 赫兹时的情况。
. use "https://www.stata-press.com/data/r17/audiometric", clear // 载入数据
. sum lft* rght* // 观察统计性质
. correlate lft* rght* // 相关系数矩阵
. pca lft* rght* // 对所有变量进行 pca 降维
Principal components/correlation Number of obs = 100
Number of comp. = 8
Trace = 8
Rotation: (unrotated = principal) Rho = 1.0000
--------------------------------------------------------------------------
Component | Eigenvalue Difference Proportion Cumulative
-------------+------------------------------------------------------------
Comp1 | 3.92901 2.31068 0.4911 0.4911
Comp2 | 1.61832 .642997 0.2023 0.6934
Comp3 | .975325 .508543 0.1219 0.8153
Comp4 | .466782 .126692 0.0583 0.8737
Comp5 | .34009 .0241988 0.0425 0.9162
Comp6 | .315891 .11578 0.0395 0.9557
Comp7 | .200111 .0456375 0.0250 0.9807
Comp8 | .154474 . 0.0193 1.0000
--------------------------------------------------------------------------
上表是 Stata 给出的解释方差表，第二列给出了每个主成分对应的特征值。该特征值的大小表示了该主成分解释数据集方差 (信息量) 的大小，所有特征值加起来就是原数据集总的方差 (信息量) 大小。该数据集总方差为 8，第一个主成分方差约为 3.93，解释了总方差的 49%(3.93/8)。同理，第二个主成分解释了总方差的 20%(1.6/8)。
每个主成分解释的百分比对应于解释方差表的第四列。第五列第 n 行数据表示前 n 个主成分累计解释方差的占比。解释方差表可以帮助我们判断提取主成分的个数，目前主要有 3 种方法可以帮助大家判断提取主成分的数量，分别是：
特征值大于 1：一般来说，如果某一项主成分的特征值小于 1，那么我们就认为该主成分对数据方差的解释程度比单个变量小，应该剔除。不过缺点是如果研究结果中某些主成分的特征值十分接近 1，那么该方法对提取主成分数量的提示作用将变得不明显。对应到本文中，我们似乎应该只保留前两个主成分，但第三个主成分是否应该保留则需要其他方法进行辅助判断； 解释数据变异的比例：既往研究认为提取的主成分至少应该解释 5-10% 的数据变异或者提取的主成分应至少累计解释 70-80% 的数据方差。根据这一指标，我们认为应该提取前 3 位主成分 (前 3 位主成分累计解释 81.53% 的数据方差)。这种判断方法的不足在于比较主观，我们既可以提取 70%，也可以提取 80%，而这 10% 的比例差异往往导致提取主成分数量的不同； 陡坡图检验：陡坡图是根据各主成分对数据方差的解释程度绘制的图。下图中横轴对应主成分解释方差从大到小点的序号，纵轴对应特征值。我们通过陡坡趋于平缓的位置判断提取主成分的数量。在本研究中，第四主成分之后的数据趋于平缓，因此我们认为可以提取前四位主成分。
. screeplot // 画出陡坡图
. pca lft* rght* // 对所有变量进行 pca 降维
Principal components (eigenvectors)
------------------------------------------------------------------------------------------------------------
Variable | Comp1 Comp2 Comp3 Comp4 Comp5 Comp6 Comp7 Comp8 | Unexplained
-------------+--------------------------------------------------------------------------------+-------------
lft500 | 0.4011 -0.3170 0.1582 -0.3278 0.0231 0.4459 0.3293 -0.5463 | 0
lft1000 | 0.4210 -0.2255 -0.0520 -0.4816 -0.3792 -0.0675 -0.0331 0.6227 | 0
lft2000 | 0.3664 0.2386 -0.4703 -0.2824 0.4392 -0.0638 -0.5255 -0.1863 | 0
lft4000 | 0.2809 0.4742 0.4295 -0.1611 0.3503 -0.4169 0.4269 0.0839 | 0
rght500 | 0.3433 -0.3860 0.2593 0.4876 0.4975 0.1948 -0.1594 0.3425 | 0
rght1000 | 0.4114 -0.2318 -0.0289 0.3723 -0.3513 -0.6136 -0.0837 -0.3614 | 0
rght2000 | 0.3115 0.3171 -0.5629 0.3914 -0.1108 0.2650 0.4778 0.1466 | 0
rght4000 | 0.2542 0.5135 0.4262 0.1591 -0.3960 0.3660 -0.4139 -0.0508 | 0
------------------------------------------------------------------------------------------------------------
上表被称为旋转成分矩阵，其实就是我们上面说到的特征向量矩阵 ，因为这个矩阵让原数据集变幻了坐标轴 (也可以说是旋转了数据)，因此有了旋转成分矩阵的名字。这个表传递了这样的信息，每一列是对应主成分的特征向量。同时因为原矩阵通过右乘该成分矩阵得到了新数据，因此这张表也说明了数据变换后各主成分与原数据各变量的关系。例如，第一列表示了：
我们前面提到过，进行 PCA 之前要进行检验，KMO 检验是一个不错的选择，我们来看一下怎么解读：
. estat kmo // 进行 KMO 检验
Kaiser-Meyer-Olkin measure of sampling adequacy
-----------------------
Variable | kmo
-------------+---------
lft500 | 0.7701
lft1000 | 0.7767
lft2000 | 0.7242
lft4000 | 0.6449
rght500 | 0.7562
rght1000 | 0.8168
rght2000 | 0.6673
rght4000 | 0.6214
-------------+---------
Overall | 0.7328
-----------------------
KMO 检验主要用于主成分提取的数据情况。一般来说，KMO 检验系数分布在 0 到 1 之间。如果系数值大于 0.6，则认为样本符合数据结构合理的要求，但大于 0.8 会更好。对于单个变量的分析结果，如果系数大于 0.5，则认为单个变量满足要求；如果系数大于0.8，则认为单个变量结果很好。在本研究中，任一变量的 KMO 检验结果均大于 0.6，总体检验结果大于 0.7，结果不算太好，但仍可以接受。Kaiser (1974) 给出了一个通用的判断 KMO 检验值的法则：
在了解完上面所有信息后，我们可以生成降维后的数据了。我们选择提取前三个主成分，通过原数据生成三个主成分 f_1、f_2、f_3。
. pca lft* rght*, components(3) // 选择前三个主成分降维
. predict f1-f3 // 生成新数据
至此我们已经完成了 PCA 的全部过程。在进行实证研究或机器学习中，经过 PCA 降维后的数据往往会比原数据在模型上拟合精度更高，也更有利于在繁杂的数据中快速找到规律，这其实就是一个数据净化的过程。需要注意的是，PCA 是需要提前对数据进行中心化的，程序里的函数已经自动执行这一步骤，各位读者如果选择自行处理的话要记得加上数据中心化的步骤。
理论 + 实证：从「读懂模型」到「折腾模型」
🎦 理论模型构建专题
📅 2022 年 4 月 23-24 日 (周六-周日)
🔑 郭凯明副教授 (中山大学)
🍓 课程主页：https://gitee.com/lianxh/emodel
4. 参考文献
Géron A. Hands-on machine learning with Scikit-Learn, Keras, and TensorFlow: Concepts, tools, and techniques to build intelligent systems[M]. " O'Reilly Media, Inc.", 2019. -PDF- Jackson J E. A user's guide to principal components[M]. John Wiley & Sons, 2005. -PDF- Abdi H, Williams L J. Principal component analysis[J]. Wiley interdisciplinary reviews: computational statistics, 2010, 2(4): 433-459. -PDF- Kaiser H F. An index of factorial simplicity[J]. psychometrika, 1974, 39(1): 31-36. -PDF- Lamberson R H, McKelvey R, Noon B R, et al. A dynamic analysis of northern spotted owl viability in a fragmented forest landscape[J]. Conservation Biology, 1992, 6(4): 505-512. -PDF- Strang G. Linear algebra and its applications[M]. Belmont, CA: Thomson, Brooks/Cole, 2006. -PDF- Mitchell M N. A visual guide to Stata graphics[M]. Stata Press, 2008. -PDF- Pearson K. LIII. On lines and planes of closest fit to systems of points in space[J]. The London, Edinburgh, and Dublin philosophical magazine and journal of science, 1901, 2(11): 559-572. -PDF- Strang G, Strang G, Strang G, et al. Introduction to linear algebra[M]. Wellesley, MA: Wellesley-Cambridge Press, 1993. -PDF- QR algorithm -Link- SPSS超详细教程：主成分分析-医小咖的文章 -Link- 再谈协方差矩阵之主成分分析 -Link- 主成分分析原理详解 -Link-
5. 相关推文
Note：产生如下推文列表的 Stata 命令为：
lianxh 主成分, m
安装最新版lianxh命令：
ssc install lianxh, replace
专题：计量专题 主成分分析-交互固定效应基础：协方差矩阵的几何意义 专题：Stata绘图 Stata可视化：biplot一图看尽方差、相关性和主成分 专题：倍分法DID Stata：平行趋势不满足？主成分DID来帮你！- pcdid
New！ Stata 搜索神器：
lianxh和songblGIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all
🍏 关于我们
连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。 直通车： 👉【**百度一下：**连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。